options(scipen=20)

library(MASS)
library(WhatIf)
library(RColorBrewer)
library(simcf)
library(verification)
library(foreign)
library(nnet)
library(xtable)
library(apsrtable)
library(ggplot2)
library(dplyr)
library(stargazer)
library(haven)

col <- brewer.pal(3,"Dark2")

data <- read_dta("cas_chimet.dta")

# Used to make Figure A1 in the Appendix

install.packages('interflex', type = "source", repos = 'http://cran.us.r-project.org') 
library(interflex)
attach(data)
# diagnostic tests for linearity (not applicable) and common support

inter.raw(Y = "act_index", D = "custodial", X ="org", data=data,
          "Ylabel"="Participation",Dlabel="Custodial",
          Xlabel="CBO Contact")

inter.raw(Y = "act_index", D = "questioned", X ="org", data=data,
          "Ylabel"="Participation",Dlabel="Custodial",
          Xlabel="CBO Contact")

data1<-data[ which(data$custodial==1), ]
ap1<-ggplot(data1, aes(x=org, y=act_index))+
  geom_jitter(shape=1,width = 0.08, height = 0.1)+
  geom_smooth(method = 'loess')+
  labs(y="\nNonvoting Participation\n",x="CSO Connection",
       title="Among those under correctional control")

ggsave(ap1, filename = paste("commonsupport_cust", ".pdf", sep = ""),width = 7.5, height = 5)


data2<-data[ which(data$questioned==1), ]
ap2<-ggplot(data2, aes(x=org, y=act_index))+
  geom_jitter(shape=1,width = 0.08, height = 0.1)+
  geom_smooth(method = 'loess') +
  labs(y="\nNonvoting Participation\n",x="CSO Connection",
       title="Among those who have been detained")

ggsave(ap2, filename = paste("commonsupport_que", ".pdf", sep = ""),width = 7.5, height = 5)

# Models hold across multiple specs (Table A9-A11 of the Appendix)
# poisson
library(AER)
m1<-glm(act_index~custodial*org+questioned*org, 
        family="poisson",data=data,weights=hweight)

dispersiontest(m1)

m2<-glm(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1, 
        family="poisson",data=data,weights=hweight)

dispersiontest(m2)
m3<-glm(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1+
          polint+eff2+dem+ind+married+unemp, 
        family="poisson",data=data,weights=hweight)

dispersiontest(m3)
stargazer(m1,m2,m3, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Correctional Control","CSO Connection",
                             "Detained by Police",
                             "Black","Latino","Other Race",
                             "Education","Female","Age: 18-34","Age: 65+","Income",
                             "Political Interest","Political Efficacy",
                             "Democrat","Independent","Married",
                             "Unemployed","Correct*CSO","Detained*CSO",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.1,0.05,0.01),
          column.labels = c("Model 1","Model 2","Model 3"),
          dep.var.labels   = "")

# regression

m1<-lm(act_index~custodial*org+questioned*org, 
        data=data,weights=hweight)

m2<-lm(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1, 
        data=data,weights=hweight)

m3<-lm(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1+
          polint+eff2+dem+ind+married+unemp, 
        data=data,weights=hweight)
summary(m3)
stargazer(m1,m2,m3, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Correctional Control","Detained by Police",
                             "CSO Connection","Black","Latino","Other Race",
                             "Education","Female","Age: 18-34","Age: 65+","Income",
                             "Political Interest","Political Efficacy",
                             "Democrat","Independent","Married",
                             "Unemployed","Correct*CSO","Detained*CSO",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.1,0.05,0.01),
          column.labels = c("Model 1","Model 2","Model 3"),
          dep.var.labels   = "")

# negative binomial doesn't converge. Quasi-poisson instead

m1<-glm.nb(act_index~custodial*org+questioned*org,data=data,weights=weight2)
m1<-glm(act_index~custodial*org+questioned*org,family=negative.binomial(theta=m1$theta),data=data,weights=hweight)

m2<-glm.nb(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1, 
      data=data,weights=hweight)
m2<-glm(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1,
        family=negative.binomial(theta=m2$theta),data=data,weights=hweight)

m3<-glm.nb(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1+
           polint+eff2+dem+ind+married+unemp,data=data,weights=hweight)
m3<-glm(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1+
             polint+eff2+dem+ind+married+unemp,family=negative.binomial(theta=m3$theta),
        data=data,weights=hweight)


m1<-glm(act_index~custodial*org+questioned*org,family=quasipoisson,data=data,weights=hweight)

m2<-glm(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1,
        family=quasipoisson,data=data,weights=hweight)

m3<-glm(act_index~custodial*org+questioned*org+black+latino+other+educ+fem+young+old+inc1+
          polint+eff2+dem+ind+married+unemp,family=quasipoisson,
        data=data,weights=hweight)

stargazer(m1,m2,m3, no.space=TRUE,type="text",
          single.row = TRUE,
          covariate.labels=c("Correctional Control","CSO Connection",
                             "Detained by Police",
                             "Black","Latino","Other Race",
                             "Education","Female","Age: 18-34","Age: 65+","Income",
                             "Political Interest","Political Efficacy",
                             "Democrat","Independent","Married",
                             "Unemployed","Correct*CSO","Detained*CSO",
                             "Constant"),
          dep.var.labels.include = FALSE,
          model.numbers=FALSE,
          star.char=c("*","**","***"),
          star.cutoffs = c(0.1,0.05,0.01),
          column.labels = c("Model 1","Model 2","Model 3"),
          dep.var.labels   = "")






